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ABSTRACT 

A technique is described that is used to improve the detection of radio-frequency interference in astronomical radio observatories. It 
is applied on a two-dimensional interference mask after regular detection in the time-frequency domain with existing techniques. The 
scale-invariant rank (SIR) operator is defined, which is a one-dimensional mathematical morphology technique that can be used to 
find adjacent intervals in the time or frequency domain that are likely to be affected by RFI. The technique might also be applicable in 
other areas in which morphological scale-invariant behaviour is desired, such as source detection. A new algorithm is described, that 
is shown to perform quite well, has linear time complexity and is fast enough to be applied in modem high resolution observatories. 
It is used in the default pipeline of the LOFAR observatory. 
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1. Introduction 

Modern telescopes in radio astronomy, such as the Expanded 
Very Large Array (EVLA) and the Low-Frequency Array 
(LOFAR), are sensitive devices that observe the sky with enor- 
mous depth and detail. The observed bandwidth of telescopes 
has dramatically increased over the last decades, and often over- 
laps with parts of the radio spectrum that have not been reserved 
for radio astronomy. Simultaneously, the radio spectrum is be- 
coming more crowded because of technological advancement. 
Therefore, radio observations are affected by man-made radio 
transmitters, which can be several orders of magnitude stronger 
than the weak celestial signals of interest. This kind of inter- 
ference, which seriously disturbs radio observations, is called 
radio-frequency interference (RFI). 

Numerous techniques have been suggested to perform the 
challenging task of RFI mitigation. They include using spatial 
information to null directions, provided in interferometers or 
multi-feed systems (Ellingson & Hampson 2002 , |Leshem et aT] 



2000; |Smolders & Hampson||2002| |Kocz et al.||2010); remov 
ing the RFI by using reference antennas (Barnbaum & Bradley 



1998:); and blanking out unlikely high values at high time reso- 
lutions (Baan et al. 2004; Leshem et al. 2000; Niamsuwan et al. 



2005; [Weber et aL|1997) . Despite the numerous possible tech- 
niques, almost any observation needs to be post-processed due 
to RFI effects. The most used technique for such a final pro- 
cessing step consists of detecting the RFI in time, frequency 
and antenna space, and ignoring the contaminated data in fur- 
ther data processing. This step is often referred to as "data flag- 
ging". Historically, this step was performed by the astronomer. 
However, because of the major increase in resolution and band- 
width of observatories, leading to observations of tens of ter- 
abytes, this is no longer feasible. The tendency is therefore to im- 
plement automated RFI flagging pipelines in the observatory's 



pipeline. Examples of these are the RFI mitigation pipeline used 
for the Effelsberg Bonn HI Survey ( Floer et al.||2010| l and the 
AOFlagger pipeline (Offringa et al.|2010a] [ 



1.1. RFI detection 

RFI detection often involves thresholding based on amplitude 
( |Winkel et al.|2006[pffringa et al.|2010b| l, although also higher 
order statistics such as the kurtosis have been used ( |Gary et ah] 
2010). The latter requires storing both the mean powers and the 
squared powers, thereby doubling the data rate, and hence is not 
always usable. Most interfering sources radiate either in a con- 
stant small frequency range, or produce a broadband peak in a 
short time range. Examples of such interferences are respectively 
air traffic communication and lightning. Consequently, an inter- 
fering source tends to affect multiple neighbouring samples in 
the time-frequency domain. These samples form straight lines, 
parallel to the time and frequency axes. An example is given 
in Fig. 1 a) which shows data from the Westerbork Synthesis 



Radio Telescope (WSRT). This line-shaped behaviour of RFI 
can be used to improve the accuracy of detection algorithms. 
An algorithm that uses this information is the SumThreshold 
method, which shows a very high detection accuracy compared 
to other methods (Offri nga et al.pOlOb) . This method is used 
in one of the steps in the AOFlagger pipeline (Offringa et al. 
2010a). An important consideration for succesful application of 
automated feature detection algorithms such as these, is that the 
signal of interest should not contain significant line-shaped fea- 
tures, as is the case with spectral line observations. Also, meth- 
ods that assume straight, one-dimensional features in the time- 
frequency domain, might not work well in situations where the 
features are curved. This can occur when both the frequency and 
the time resolutions are high enough to resolve frequency varia- 
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Fig. 1: Typical spectral line RPI received in a short period of WSRT data around 117 MHz. It is likely that such RFI sources 
transmit continuously within a small bandwidth. Panel |(a)| shows the original observation, while panel |(b)| shows what the AOFlagger 
with default settings would flag without morphology-based flagging. Detection is quite accurate, but some of the detected lines in 
panel [(b)] are not continuous. It is likely that those RFI sources were active in the gaps as well. Morphology-based detection will 
help in such cases. The plot shows Stokes I amplitudes of the cross-correlation of antennas RTO x RT1, which is a 144m East-West 
baseline. A single pixel is 10 seconds x 10 kHz of data. 



tion in sources, for example when sources are Doppler shifted 
or vary intrinsically, such as with certain radar signals. With 
LOFAR, we see very few such sources. 

Typically, the received power of interfering sources varies 
over time and frequency. This happens because of several effects, 
such as intrinsic variation of the source; changing ionosphere; 
and because of instrumental effects. A typical example of the lat- 
ter, which is present in almost every observation, is the change of 
the telescope's gain towards a terrestrial source as the telescope 
tracks a field in the sky. Like time variation, frequency variation 
can be caused intrinsically by the source. The instrument also 
adds frequency-dependent gain, for example due to imperfect 
band-pass filters. Even though a source might be continuously 
received by the telescope, thresholding detection methods might 
fail to detect the interferer over its full range due to the varia- 



tion in received power. Figure 1 b) shows an example where this 
is likely the case. Increasing the sensitivity of the thresholding 
method might help somewhat, but will also cause an increase of 
false positives. While some falsely detected samples are tolera- 
ble, they should be kept minimal in order to avoid data bias and 
insufficient resolution. 

Using mathematical morphology for RFI detection is not a 
new idea; a dilation is often used during RFI processing to flag 
areas near high values in the time-frequency domain. An exam- 
ple of this can be found in lWinkeF et al. (2006), where windows 
of 5 time steps x 5 frequency channels around detected samples 
are flagged. However, standard morphological techniques are not 
scale invariant. An operator is called scale invariant if scaling 
its input results in the same scaling of its output. An ordinary 
dilation will cause sharp RFI features to create a high amount 
of false positives, while flagging smooth RFI features requires 
a very large dilation kernel. Another scale-dependent technique 
used for RFI detection is to consider the statistics of time steps 
and frequency channels. In this paper, we will show that scale 



invariance is a desirable property of RFI detection algorithms. 
This paper provides: 

- A detailed description of a recently-introduced morphologi- 
cal technique for RFI detection; 

- Analysis of the technique and a comparison with an ordinary 
dilation, using simulations and real data from two different 
radio-observatories; 

- A novel fast algorithm with linear time complexity to imple- 
ment the technique. 

1.2. Outline 

We discuss a morphological technique that can be used to im- 
prove RFI detection. The method flags additional samples that 
are likely to be contaminated with RFI, based on the morphology 
of the flag mask output of a thresholding stage in the pipeline. 
In sect. [2] we describe the technique and show a fast algorithm 
to implement it. We present some results of the method on sim- 
ulated data and real data in Sect. [3] Finally, we summarize and 
discuss the results in Sect. [4] 



2. The scale-invariant rank operator 



RFI features such as in Fig. 2|[a)| are common in radio obser- 
vations, and can occur at different scales. However, a morpho- 
logical dilation is not scale invariant, and will thus necessar- 
ily work better for some RFI features than others. To overcome 
this problem, we will describe and analyse a morphological rank 
operator that is scale invariant Scale invariance is a desirable 
property of RFI detection algorithms, because (a) it implies the 



1 The mathematical properties of this technique will be analysed in 
more detail in van de Gronde et al., 2012, in preparation. 
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method can be applied on data with different resolutions with- 
out changing parameters and (b) the time and frequency scale of 
RFI itself can be arbitrary, so any method to detect RFI should 
work equally well for RFI at different scales. In practice, RFI 
seems to behave in a more or less scale-invariant manner at the 
resolution of LOFAR, as for example can be seen in Fig. [TJ so 
we should also use a scale-invariant method to detect it. This 
scale-invariant behaviour of RFI breaks down at high time and 
frequency resolutions, at which many features become diagonal 
in the time-frequency plane. 



Horizontal Vertical 



The proposed technique was first mentioned in Offringa et al. 
(2010a), as it is part of the AOFlagger, which is the default 



LOFAR RFI detection pipeline. In that article the operation was 
referred to as a dilation, however, it does not strictly adhere to 
all the properties of a morphological dilation. For example, we 
will see that the operator p is not distributive over the union set 
operator: p(X U 7) + p(X) U p(Y) for some X and Y. Because a 
rank operator flags points for which the number of flagged points 



in a neighbourhood exceeds a threshold ( Goutsias & Heijmans 
|2000| §3.4, |Sollle 2002), we will refer to the operator p as the 



scale-invariant rank (SIR) operator. 

In this paper, we will describe the method in-depth and anal- 



yse its effectiveness. In Offringa et al. (2010a I, it was mentioned 
that the full algorithm has a time complexity of 0(N 2 ), N being 
the input size of the SIR operator, but by making the algorithm 
less accurate, an implementation of O(NxlogN) was mentioned 
to be possible. Here, we will introduce a faster algorithm with 
linear time complexity, which is also an exact implementation 
of the SIR operator. 

2.1. Description 

Consider F, a set of positions in the time-frequency domain, 
such that a sample at time t and frequency v has been flagged 
when (t, v) € F. Assume F is the result of a statistical detection 
algorithm, such as the SuraThreshold algorithm. We will apply 
the SIR operator in time and frequency directions separately, and 
define the sets 0, and cD v to contain the flags of a slice in time 
and frequency direction: 



0, = {(i,v)eF| s = t}, 
Q> y ={(t,fi)€F\n = v}. 



(1) 
(2) 



A single one-dimensional set 0, or cD v is the input for the SIR 
operator. The operator considers a sample to be contaminated 
with RFI when the sample is in a subsequence of mostly flagged 
samples. To be more precise, it will flag a subsequence when 
more than (1 - rj)N of its samples are flagged, with N the number 
of samples in the subsequence and r\ a constant, < r\ < 1 . Using 
p to denote the operator, the output p(X) can be formally defined 



as 



p(X) = (J {[71, 72) 



(3) 



#(Xn [71,72)) > (1 -77)(72- 71) , 



with [7i, Y2) a half-open interval of 0, or (D v , and the hash sym- 
bol # denoting the count-operator that returns the number of el- 
ements in the set. In words, Equation [3 defines p(X) to consist 
of all the samples that are in an interval [71, 72), in which the 
ratio of samples in the input X is greater or equal than (1 - 77). 
Parameter 77 represents the aggressiveness of the method: with 
77 = 0, no additional samples are flagged and p(X) = X. On the 
other hand, 77 — 1 implies all samples will be flagged. Figure [2] 
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(a) Input (b) Union 




— — 



(c) Horizontal (d) Vertical first 
first 



Fig. 3: Example outputs of the SIR operator in which the one- 
dimensional output has been combined in three different ways. 
Panel (a) is the input, panel (b) shows the result of performing a 



union on the outputs of both directions, and in panels (c) and (d) 



the SIR operator was first applied in, respectively, the horizonta 
and vertical direction. Parameter 77 was 0.5 in this example. 



shows an example of a simulated Gaussian broadband RFI fea- 
ture, and the input and output of the SIR-operator. 

The one-dimensional outputs can be remapped to the original 
two-dimensional domain in various ways. A simple and useful 
way is to perform a logical union of ©J = p(0,) and <&' v - p(<J> v ), 
the flags on respectively the time and frequency outputs: 



F' 



(4) 



An alternative is to initially apply the SIR operator only in 
one direction, i.e., on the sets that correspond with either the 
time or frequency direction, and subsequently applying the SIR 
operator on the outputs of the first in the other direction. The lat- 
ter is more aggressive than the former. The result also depends 
on which direction is processed first. The difference is demon- 
strated in Fig. [3] and an example of how that would work out 
on actual data is given in Fig. |4] Optionally, the operator can be 
applied in frequency and time directions with different 77, if one 
suspects that RFI acts differently in either direction. 



2.2. Properties & parameters 

Consider the case in Equation Q when a subsequence of ar- 
bitrary length is flagged. Since the fraction of flagged samples 
within the subsequence is explicitly used to define its output, the 
operator is scale invariant. Formally, an operator p is scale in- 
variant if and only if p{AX) = Ap{X), i.e., scaling the input X 
with A followed by p is equal to scaling the output p(X) with A. 
We will now give a formal proof of the scale invariance of the 
SIR operator. 

Proof. With p the SIR operator, we will scale the input X with 
factor A. If A - Owe trivially have that p{AX) = Ap(X). Also, if 
p(AX) = Ap(X) for A > 0, it is not difficult to see that we also 
have p(-AX) = -Ap(X), as mirroring the input will mirror the 
output. Therefore, assume without loss of generality that A > 0. 
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Fig. 2: Simulation of a typical broadband RFI feature with Gaussian frequency profile as used in the ROC analysis. Panel (a): isolated 
RFI feature; panel (b): when noise is added, a part of the feature becomes undetectable; panel (c): flagged with the SumThreshold 
method; panel (d): with SIR operator applied, parameter rj = 0.2. 




(a) Original data (b) Intersection 




(c) Union (d) Time first 




(e) Frequency first (f) Union of (d) and (e) 



Fig. 4: Example of the SIR operator applied on a LOFAR observation, displaying five different methods to make the SIR operator 
two-dimensional. The visibilities shown are from baseline CS003 x CS007 of a LOFAR low-band-antenna (LBA) observation with 
3s x 0.8 kHz resolution. This observation part was selected as an example because it has a two-dimensional RFI structure. Such 
RFI is less common, hence this is not a typical case. With the exception of the intersection, there is no difference between the 
different methods on the thin lines below 32.1 MHz. Applying the operator sequentially (panels d, e and f) is more aggressive for 
the two-dimensional structures, as it will flag samples that have diagonal neighbours that are flagged. Intersecting the two methods 
(panel b) will only flag concave samples. Pink is pre-flagged by the SumThreshold method, yellow is added by the SIR operator. 
A value of 77 = 0.2 was used in this example. 
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Now, substituting X with AX in Equation ([3]) results in 

p(AX) = |J {[71, 72) | 

#(AX n [71, 72)) > (1 - T])(Y2 - 71)} . 

By using Z\ -Y\/A and Z2-Y2/A, this can be rewritten to 

p(AX) = [J{[AZ1,AZ2) I 

#(,IX n [AZl,AZ2)) > (1 - t?)(/IZ2 - ^Zl)} . 

If we assume continuous positions, both the left side and the 
right side of the comparison can be scaled by l/A: 

p(AX) = yj{[AZl,AZ2) I 

f(ln [Z1,Z2)) > (1 - j])(Z2 - Zl)} , 

and by using [AZ\,AZ2) - A[Z\,Z2) and the definition in 
Equation ([3j, this is equivalent to p(AX) = Ap(X). 

Because the time and frequency dimensions are obviously dis- 
crete and finite when applied on radio observations, in practice 
the scale invariance is limited by the resolution and size of the 
data. 

The aggressiveness of the SIR operator can be controlled 
with the 77 parameter, which can be chosen differently for the 
time and frequency directions. Because the method is scale in- 
variant, the choice of r\ can be made independent of the time and 
frequency resolutions of the input. The default 77 parameter in 
the LOFAR pipeline is currently 77 = 0.2 and is equal in both 
directions. This value has been determined by tweaking of the 
parameter and data inspection, e.g. by looking at the resulting 
time-frequency diagram and projections of the data variances. 
The results were checked for many observations. Higher values 
seem to remove too much data without much benefit, while some 
RFI is left undetected with lower values. The value works well 
for various telescopes and on different time and frequency reso- 



lutions. We will evaluate this setting in Section 3.2 



Since most telescopes observe two linear or circular polar- 
izations, RFI detection can consider each (correlated) polariza- 
tion individually, and the operator can be applied on each pro- 
duced mask independently. However, the flag masks are often 
kept equal between the polarizations, because calibration might 
become unstable when, for a particular sample, part of the polar- 
izations are missing. Moreover, if one of the polarization feeds of 
the telescope has been affected by RFI, it is likely that the others 
also have been affected. For these reasons, the approach taken in 
the LOFAR pipeline is to use the SumThreshold method on all 
four cross-correlated polarizations (XX, XY, YX and YY) indi- 
vidually, then flag any sample for which at least one polarization 
has been flagged, and finally apply the SIR operator once on the 
combined mask. 



2.3. The algorithm 

A straightforward implementation of the operator in 
Equation <[3j is to test each possible contiguous subsequence. 
In this case, if N is the number of samples in the sequence 0, 
or ct> v , 0(N 2 ) sums of subsequences have to be tested. Since 
the sums of all subsets can also be constructed in quadratic 
time complexity, the total time complexity of a straightforward 
implementation is 0(N 2 ). We will now show an algorithm that 
solves the problem with linear time complexity. The algorithm 
is somewhat similar to the maximum contiguous subsequence 
sum algorithm. 



Listing 1: Linear time complexity algorithm for the scale-invariant 
rank operator 

function ScalelnvariantRankOperator 
Input : 

N : Size 

Q : Input array of size A' 

(fl[i] = 1 => i is flagged, 
Q[i] = otherwise) 
: Aggressiveness parameter 



Output : 



Output flag array of size N 



1 : begin 

// Initialize *F 
2: for x = 0...N- 1 do ^[x] 



1 + nw 



18 

11 
12 

13 
14 



15 
16 
17 
18 
19 
2Q 
21 



// Construct an array M such that: 
// M{x) = 2 je{0...x-l}: <9[j] 
M[0]<- 

for x = 0...N- 1 do M[x + 1] «- M[x] + ?W 

// Construct array P such that: 
// M[P[x]] = min M [j] : 0< j<x 
P[0]<- 

for x = 1 . . . N - 1 do 
P[x]<-P[x-1] 

if M[P[*]] > MM then P [*]<-* 
end for 

// Construct array Q such that: 
// M[Q[x]] = max M [j] : x< j<N 
Q[N- l]<-iV 
for x = N- 2...0 do 
fiM«-GDc + l] 

if MlQ\_x~\~\ < Mix +12 then x+ 1 

end for 

// Flag sample x if M[Q[x]] - M[P[x]] > 
for x = . . . N - 1 do 

if M[g[x]]-M[P[x]]> then 

n' [x] «- 1 

else 

n' [x] «- 

end if 
end for 



22: return fi' ; 
23: end 

Listing 1 shows a direct algorithm to solve the SIR operator 
problem. 

Proof. Using the definition of D.(x) and Q.'(x), such that 1 in- 
dicates that x is flagged and that it is not, we can rewrite 
Equation Q as 



Q'(x) 



1 if 37i < x, 7 2 > x, such that 

y=Yi 

otherwise. 



(5) 



In line 2, the array ^(y) is initialized such that Tfjy) = 77 in case 
y is flagged, and *P(y) = 77-1 otherwise. Equation <|5j can now 
be rewritten to the following test: 



Q'(x) 



1 if 37 <x, 37 2 > x : F £ V(y) > 

y=Y, (6) 

otherwise. 
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Line 3-4 initialize M(x) for < x < N to 

x-i 



7=0 



so that Equation |6]) can be rewritten as 



Q'(jc) = 







if 3Yi <x,3Y 2 >x: 

M(Y 2 )-M(Y l ) > 
otherwise. 



(7) 



Because we are only interested in £2'(x) in the range < x < N, 
we can limit the search for Y\ and Y 2 to < Y\ < x < Y 2 < N. 
There exists Y\ and Y 2 in this range such that M{Y 2 )-M(Y\) > 0, 
if and only if 



max M(y) - min M(y) > 0. 

v:.v<v<JV v:0<v<.v 



Lines 5-14 make sure that P and Q are initialized for < x < N, 
such that 



P(x) = argmin M(y), 

y£0...x 

Q(x) = argmaxM(y). 

y©c+l.,JV 



Finally, this allows Equation |7]) to be rewritten as 



Q!{x)= \ l ifM(Q(x))-M(P(x))>0 
1 otherwise, 



which is performed and returned in lines 15-23. 

The algorithm is &(N), and performs 3iV additions or sub- 
tractions and 3N - 2 comparisons on floating point numbers. 
The algorithm uses the temporary arrays W, M, P and Q, each of 
size N, with the exception of M which is of size N + 1 . Array *F 
can be optimized away and the input O can be reused for output 
by assigning directly to it in lines 17 and 19. The total amount of 
temporary storage required is thus about YV floating point values 
and 2N index values, thus O(N). When the function is applied 
on a two-dimensional image, as in the case of RFI detection, 
the temporary storage is negligible, as the number of processed 
slices is usually much larger than one or two. If rj is expressed as 
a ratio of two integer values, it is possible to scale all values and 
only use integer math. 

The algorithm has been implemented in C++ and takes 
around 40 lines of cod^f] 

Because the problem is somewhat similar to the maximum 
contiguous subsequence sum (jBentley 1984) and the all max- 
imal contiguous subsequence sum problems, it might be pos- 
sible to parallelize the algorithm by similar means, e.g. as in 
( jAlves et al.||2005| ). Moreover, parallel algorithms exist for the 
prefix sum/min/max calculations. For the specific application of 
RFI detection for LOFAR, the pipeline has already been max- 
imally parallelized by flagging different baselines and/or sub- 
bands concurrently. Unlike parallelizing on the algorithm level, 
this requires no communication between the different processes. 



2 The implementation is part of the AOFlagger and can be down- 
loaded from http : //www . astro . rug . nl/rfi- software 
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Fig. 5: Computation time versus input size with the different al- 
gorithms and fixed 77 = 0.2. The average over 1000 runs was 
taken for each different configuration. 



3. Analysis & results 

3.1. Performance 

Figure [5] displays the performance of the C++ implementa- 
tions and compares the linear algorithm with the approximate 
0(N log N) algorithm and the full quadratic algorithm. The mea- 
surements have been performed on a regular desktop with a 3.07 
GHz Intel Core i7 CPU, using only one of its cores. The time 
complexities of the three algorithms for increasing N behave 
as expected. The linear algorithm is faster in all cases, even for 
small N. The 0(N log N) time complexity algorithm is more than 
one order of magnitude slower at both small and large N. The lin- 
ear algorithm has been executed with different values for rj and 
the results are shown in Fig. [5] Except for some slight variations 
— especially for 77 = — the algorithm's speed is independent 
of 77. 

In the LOFAR pipeline, it takes 3.8 seconds to process a 
single sub-band for a single baseline, assuming 100,000 time 
steps and 256 channels (which is common). Of these 3.8 sec- 
onds, only 49 milliseconds (1.3%) are spent applying the SIR 
operator. In common applications, an observation contains on 
the order of a 1,000 baselines and 250 sub-bands. The pipeline 
is heavily parallelized by concurrently flagging baselines over 
multiple cores and sub-bands over multiple cluster nodes. In this 
case, the pipeline's performance is dominated by disk access, 
and the relative contribution of the SIR operator is even smaller. 

3.2. Accuracy analysis 

The performance of the SIR operator was tested by using re- 
ceiver operating characteristics (ROC) analysis. To do so, a 
ground truth needs to be available, which can only be accu- 
rately acquired in a simulated environment. As discussed pre- 
viously, a very large fraction of RFI is line-like. The samples 
on such a line are not uniform due to intrinsic effects or instru- 
mental gain variations. Therefore, we have used simulations of 
four line-shaped RFI features as displayed in Fig. |6j (a) a sin- 
gle Gaussian that reaches its 3<x point at both borders and is 1 
in the centre; (b) three periods of a sinusoidal function which 
is scaled between zero and one; (c) the Gaussian feature, but 
slanted by 1/50 fraction; and (d) a burst-like signal in which the 



6 



A.R. Offringa et al.: A morphological algorithm for improving RFI detection 




Fig. 6: Features used for the accuracy analysis. Panel (a): Feature with Gaussian slope; panel (b): Sinusoidal feature; panel (c): 
Slanted feature with Gaussian slope; panel (d): Burst feature with samples drawn from a Rayleigh distribution. 



amplitude levels are drawn from a Rayleigh distribution with 
mode cr = 0.6. All features are three samples wide. Complex 
Gaussian distributed noise with cr - I was added to the image, 
such that the amplitudes are Rayleigh distributed. The created 
two dimensional image of size 180 x 1024 was subsequently 
flagged by the SumThreshold method with settings as in the 
LOFAR AOFlagger pipeline, and the created flag mask was used 
as input. 

To estimate the performances, the true and false positives ra- 
tios (TP and FP ratios respectively) were calculated after detec- 
tion. We created a fuzzy ground truth mask in which a value of 
one corresponds with maximal RFI, zero corresponds with sam- 
ples not contaminated by RFI, and values in between correspond 
with lower levels of RFI contamination. Fig. 1 a) shows for ex- 



ample the ground-truth mask of the Gaussian feature. Given a 
sample with ground truth value /3, if the corresponding sample 
was flagged by the method, it would be counted with ratio /3 as 
a true positive and 1 - /3 as a false positive. The total TP and FP 
ratios are the sum of all the TP and FP values, divided by the 
total sum of positives and negatives in the test set, respectively. 

The SIR operator and a standard morphological dilation 
have been applied in the direction of the feature, i.e., verti- 
cal/frequency direction. The true and false positives were varied 
by changing the parameter rj or the dilation size for respectively 
the SIR operator and the dilation. Different runs gave slightly 
different results because of the introduced Gaussian noise, hence 
the simulation was repeated 100 times and the results were aver- 
aged. 

Figure [7] shows the average results. The shadowed areas in 
panels [(b)! and |(c)| show the standard deviation over the 100 runs. 
In the case of the Gaussian RFI feature, the SumThreshold pre- 
detection removes on average 91.3% RFI power, while simulta- 
neously falsely flagging a ratio of 0.38%. Hence, if the methods 
do not flag any additional samples, they have a TP/FP ratio of 
91.3%/0.38%, and this is therefore the start of both ROC curves 
for this RFI feature. With rj = 0.48, the SIR operator flags all the 
RFI features with 100% TP with a FP ratio of 1 .36%, with the ex- 
ception of the slanted feature. The SumThreshold pre-detection, 
dilation and SIR-operator work less well on the slanted feature, 
and fail to detect it with 100% even at very high sensitivity. The 
dilation operator needs a size of 32.8% of the height of the image 
to detect the vertical RFI features. Since it will dilate any falsely 
detected input sample equally, its FP ratio is 46% with this set- 
ting. Changing the signal-to-noise ratio (SNR) of the features 



changes the scaling of the ROC curves, but the relative differ- 
ence between the two methods remains the same. 

Given the various types of RFI, Fig.[7]shows that (I) the SIR 
operator is extremely accurate on straight features, by detect- 
ing all previously undetected samples with only a very slight 
increase in false positives; (II) the SIR operator is superior to the 
dilation in all tested situations; and (III) a setting of 77 ~ 0.2-0.4 
seems to be a good compromise between FP and TP ratios. 

These tests have been performed by applying the operators 
in one dimension. When applied in two dimensions by using the 
output of the first dimension as input for the second, the compar- 
ison between the dilation and the SIR operator will diverge even 
more, because the false positives created by the first dimension 
will be multiplied by the repeated application towards the sec- 
ond dimension. An example of a two-dimensional application 
is shown in Fig. [8] Certain RFI sources create more complex 
shapes in the time-frequency domain, and contaminate larger 
non-line like areas. These RFI sources cause higher values in 
the output of the Gaussian smoothing, which is commonly part 
of the earlier RFI detection stage, and consequently some of the 
lower RFI levels of the RFI feature are not flagged. We have 
seen that the SIR operator will work very well on such features, 
because it fills the feature and slightly extends the flags in all 
directions. 

It should be noted that one of the assumptions made for the 
SIR operator to improve detection accuracy, is that parts of the 
RFI features are not detectable by amplitude thresholding. In 
practice, however, a small subset of received RFI sources does 
contribute to an observation with sufficient strength to detect the 
entire feature with amplitude thresholding. Such transmitters are 
the worst-case situation for the SIR operator, as the operator will 
enlarge the flag mask relative to its length, but any samples it 
flags extra are false positives. Note that it is not useful to per- 
form ROC analysis of such a situation, as the true positives will 
be constant. The number of false positives can easily be calcu- 
lated, and scales linearly with 77 and the duration of the transmit- 
ter. For example, when applying the SIR operator with r\ = 0.2 
on a strong RFI transmitter that occupies ten minutes of data in 
one channel, the operator will falsely flag two minutes of the 
channel before and after it. An example of a band that contains 
intermittent transmitters is the air traffic communication band of 
1 18-137 MHz. Nevertheless, while some of these transmitters 
are indeed strong, e.g., when they fly through a beam sidelobe, 
there are also many transmitters at this frequency that are too 
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(left axis); dashed lines: false positives (right axis). 



Fig. 7: Analysis of the receiver operating characteristics of the SIR operator and a standard dilation on simulated data. Marks (a)-(d) 
correspond with the features shown in Fig. [6] respectively a Gaussian broadband feature, a sinusoidal feature, a slightly slanted 
Gaussian feature and a burst-like feature. The shadowed areas show lcr levels over 100 runs. 



weakly received to be detected all of the time. Consequently, 
some of them are only partially detected with amplitude thresh- 
olding. This is why we expect better results using the SIR oper- 
ator even in these bands, compared to using a dilation. 



Figure [9] shows a WSRT example that contains many differ- 
ent RFI kinds. The initial SumThreshold method detects the 
RFI quite accurately, but it leaves some parts of the last 1.5 
hours unflagged. This is solved by the SIR operator, although, 
because of the sudden start of the RFI, it falsely flags about 20 
minutes of data before the start of the RFI. The strong RFI pro- 
duced by the sporadic transmitter around 140 MHz is flagged 
by the SumThreshold method, but in this case it is likely that 
these channels have been occupied all of the time. Therefore, the 
SIR operator gives the desired result by increasing the flags in 
those channels. All in all, the baseline might be somewhat over- 
flagged. Nevertheless, it does allow further data reduction with- 
out manual intervention, and without thresholding part of the 
noise. Moreover, this case is exceptional, and the sudden start of 
very strong broadband RFI is (fortunately) seldomly seen, while 



the sporadic transmitters such as the one at 140 MHz are seen 
very often. 

A final remark on the ROC analysis performed here is that 
the given absolute true/false positive ratios are not an accurate 
representation of actual RFI detection, because our two models 
are very simplistic and based on the assumption that RFI be- 
haves in a well defined manner. Establishing absolute true/false 
positive ratios would require a detailed statistical model of the 
behaviour of RFI. A realistic estimate for the number of samples 
occupied by RFI with LOFAR is in the order of a few percent 
( |Ofiringa et al.|201 



4. Conclusions and discussion 



From panel Q a) it is clear that the SIR operator is much more 
suitable to detect the tested kinds of RFI than an ordinary mor- 
phological dilation. A value of 77 = 0.2 was determined by 
tweaking and validating the results to be a reasonable setting 
for the LOFAR RFI pipeline, and has been used in the default 
LOFAR pipeline for over a year. Panel |^[b)| shows that this 
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(a) Horizontal SIR operator 



(b) Vertical SIR operator 



(c) SIR operator in both directions 
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Fig. 8: Gray-scale plots showing examples of the effectiveness of two morphological techniques on the data from Fig. [T] The 
pink samples have been set by the SumThreshold algorithm and the yellow samples have additionally been detected with the 
morphological techniques. Panels a-c show results of the SIR operator with r\ — 0.2 in the time direction and/or r\ = 0.3 in 
frequency direction. Panels d-f show an ordinary dilation with a horizontal kernel of five pixels and/or a vertical kernel of three 
pixels. 



agrees approximately with the simulations: at 77 = 0.2, the ver- 
tical features have almost been completely detected by the SIR 
operator (Gaussian: 98.9%, a 7.6% increase, sinusoidal: 99.9%, 
5.8% increase, burst: 100%, 6.7% increase) with a minor in- 
crease in the false positive ratio (Gaussian: 0.69%, an 0.31% 
increase, sinusoidal: 0.95%, 0.42% increase, burst: 1.3%, 0.08% 
increase). The slanted feature is not as accurately detected (86%, 
6.1% increase), but the method does enhance the detection. It is 
hard to give a similar optimal value for the dilation operation, 
since the false positives scale linearly with the size of the dila- 
tion kernel. Therefore, it depends on what is an acceptable loss 
in terms of the false positives. 

In the case of simulated Gaussian broadband features, 
only 8.7% of the RFI power was not detected by the 
SumThreshold method. For the sinusoidal and burst features, 
the SumThreshold method performs even better. Therefore, the 
total benefit of the SIR operator might seem small. However, we 
think that there are strong reasons to use the method: 

- The added false positives are almost negligible, and the 
chances of biasing your data are much smaller compared 
to using amplitude thresholding exclusively. For example, 
thresholding biases the final distribution of uncorrelated 



white noise, while morphologically extending a flag mask 
does not. For these reasons, it is preferable to use morphol- 
ogy to find the final few RFI samples, compared to lowering 
the threshold. 

- The method is extremely fast and simple, and its processing 
time is almost negligible in a full RFI pipeline. 

- We have seen situations in which even the low ratio of false 
negatives that are leaked through an amplitude-based RFI de- 
tection pipeline can cause calibration to fail. Empirically, we 
have seen an improvement of the calibratability of LOFAR 
observations by using the morphological method. 



Section 3.2 describes that strong intermittent (on -minute 



scale) RFI transmitters are probably the worst case for the SIR 
operator, as in these cases the application of the operator with 
77 = 0.2 could in theory yield 40% false positives. However, 
because of LOFAR's high resolution, in combination with the 
SumThreshold's unprecedented detection accuracy, the total 
percentage of flagged data in the case of LOFAR is only a few 
percent. This implies that even if a large ratio of these were 
strong intermittent transmitters - which is unlikely - the benefits 
of not having to manually consider data quality in cases where 
the technique does help, probably outweigh the ~1% added false 
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Fig. 9: Example of an interesting but uncommon WSRT case: part of a baseline of an observation at 140 MHz that suffered unusually 
strong broadband RFI during the last 1.5 hrs. It also contains many different kinds of transmitters that mostly occupy constant 
channels. The vertical stripes are fringes of celestial sources, hence contain the information of interest. The image shown is 2000 x 
250 samples in size. 



positives. If it turns out that some bands do have mostly strong 
transient transmitters, the r] parameter could become a function 
of frequency. At the moment, application of the SIR operator 
seems to be helpful at any frequency. 

In this paper we have assumed scale-invariant behaviour for 
RFI. In reality this might not be entirely accurate, so instead 
of using a threshold that grows linearly with the scale, as in 
our definition of the SIR operator, it might be better to have a 
threshold that depends on the scale in a non-linear fashion. Also, 
when looking at the problem from a statistical point of view, RFI 
might not be equally likely to occur on all scales. For example, 
RFI might be less likely to occur on a scale of days than on a 
scale of seconds. When it does occur on large scales, it is doubt- 
ful that we actually need to extend the detected intervals in a 
scale-invariant manner, because the signal would likely already 
be detected at a smaller scales and gaps would likely be filled. 
Such considerations would suggest that it might be better to have 
a threshold that grows less than linearly for large scales. Better 
RFI statistics and RFI modelling might provide the required in- 
formation for assessing such considerations. 

Several options are available to apply the SIR operator on a 
two-dimensional input. As shown in Fig. [4] the intersection of 
the results in both directions does not extend line RFI, thus is 
not useful in this context. A union does extend such RFI, but 
does not extend the flags diagonally. Processing the directions 
sequentially might therefore be beneficial for RFI that has struc- 
ture in both frequency and time, as this kind of RFI does likely 
also slightly contribute in the diagonal direction. The difference 
between processing time first, frequency first or taking the union 
of both, is small. Taking the union overcomes the somewhat ar- 
bitrary decision of which direction to process first. In the case 
of LOFAR, we decided to only perform filtering time first, be- 
cause taking the union of both time and frequency first is more 
expensive. 

Morphology can be used in several image processing tasks, 
for example in feature detection. Often, generic morphological 
operations need to be applied on different resolutions. In such 



cases, the scale-invariant operation to extend binary masks as 
presented here might be generally useful. 

So far, we have considered combinations of one-dimensional 
application of the SIR operator in order to use it for our two- 
dimensional application in the time-frequency domain. For this 
application, but also for more generic applications, it might be 
interesting to consider a true two-dimensional version of the SIR 
operator. While the one-dimensional operator selects all sub- 
sequences (lines) with a ratio > r\ of flagged values, a two- 
dimensional operator would select all rectangles that have a ratio 
> r] of flagged values. It is however likely that such an operator 
can not be implemented with a linear time complexity, which 
makes it less attractive for the large data rate of LOFAR. 

We have shown that even slightly slanted features are harder 
to detect accurately. Fortunately, in the case of LOFAR, such fea- 
tures are very rare. If the features to be detected have a known 
orientation that is not parallel to one of the axes, it might be 
an option to apply the operator in the direction of the features. 
While a trivial implementation can apply the operator along 
fixed lines, some work might be necessary to maintain transla- 
tion invariance (Soille & Talbot 200l]l. 
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